查看原文
其他

癌基因都是肿瘤的风险因子吗

生信技能树 生信技能树 2022-08-20

在:癌基因一定在肿瘤部位高表达吗 我们探索发现并不是使用的癌基因都在肿瘤部位高表达,也不是所有的高表达基因都是癌基因,对抑癌基因也是如此。这个就很有意思的,因为癌基因的定义就是那些在肿瘤里面过度激活的基因,而抑癌基因就是在肿瘤里面失活的基因,不过过度激活不一定要转录本大量增加,可能是其它生物学机理,比如蛋白质产物大量增加,又或者说蛋白质产物效果增强。

同理,我们会问另外一个问题,就是癌基因都是肿瘤的风险因子吗,它高表达会导致癌症比如死的越来越快吗?反之,抑癌基因一定是肿瘤的保护因子吗,它表达量越高癌症病人越受到保护吗,因为想当然的我们会认为抑癌基因能抑制癌症嘛,所以它表达量越高越好。

同样的,我们可以使用TCGA数据库的多种癌症来举例说明这两个问题:

整理表达量矩阵和生存分析

首先,我们选择同样的 TCGA-CDR-SupplementalTableS1.xlsx 文件里面的生存信息,在前面的教程给出来了下载地址:

library(readxl)
phe=read_excel('../MC3/TCGA-CDR-SupplementalTableS1.xlsx',sheet = 1
survdata=as.data.frame(phe[,c(2,3,26:27)] )
survdata$OS.time =  survdata$OS.time/365
boxplot(survdata$OS.time)
survdata=survdata[survdata$OS.time>0,]
table(survdata$type)
survdata=na.omit(survdata)
table(survdata$type)

可以看到,绝大部分癌症类型里面的有生存信息的肿瘤病人数量还可以:

> table(survdata$type)

 ACC BLCA BRCA CESC CHOL COAD 
  91  410 1083  294   45  442 
DLBC ESCA  GBM HNSC KICH KIRC 
  47  185  595  527  112  535 
KIRP LAML  LGG LIHC LUAD LUSC 
 288  173  511  371  509  497 
MESO   OV PAAD PCPG PRAD READ 
  85  582  184  179  500  164 
SARC SKCM STAD TGCT THCA THYM 
 261  454  417  134  506  123 
UCEC  UCS  UVM 
 546   56   80 

接下来就是把表达量矩阵信息跟这个临床信息相对应,这里我们仍然是 使用代码批量读取前面的教程:居然有如此多种癌症(是时候开启pan-cancer数据挖掘模式),的rdata文件里面的表达量矩阵,然后根据TCGA数据库的病人的ID来对表达量矩阵信息跟这个临床信息相对应。


fs=list.files('Rdata/',
              pattern = 'htseq_counts')
fs

# 首先加载每个癌症的表达量矩阵
# 在每个癌症里面与临床信息相对应。

pre_sur  = lapply(fs, function(x){
  # x=fs[9]
  pro= gsub('TCGA-','',gsub('.htseq_counts.Rdata','',x)) 
  print(pro) 
  load(file =  file.path('Rdata/',x))  
  Group = ifelse(
    as.numeric(substring(colnames(pd_mat),14,15)) < 10,
    'tumor','normal'
  )
  table(Group)
  
  exprSet <- pd_mat[,Group=='tumor']
  dim(exprSet) 
  exprSet[1:4,1:4]
  pid=substring(colnames(exprSet),1,12)
  exprSet=exprSet[,!duplicated(colnames(exprSet))]
  colnames(exprSet) = substring(colnames(exprSet),1,12)
  
  # 1. prepare data for coxph----
  ## 批量生存分析使用coxph回归方法
  ## http://www.sthda.com/english/wiki/cox-proportional-hazards-model
  
  phe <- survdata[survdata$type==pro,c(1,4,3)]
  colnames(phe)=c('p','time''event')
  phe=phe[!duplicated(phe$p),]
  rownames(phe)=phe$p
  
  ids= intersect(  rownames(phe),   colnames(exprSet) )
  phe=phe[ids,]
  exprSet=exprSet[,ids]
 
  return( list(
    phe=phe,
    exprSet=exprSet 
  ))
  
})
names(pre_sur) = gsub('TCGA-','',gsub('.htseq_counts.Rdata','',fs)) 
do.call(rbind, 
        lapply(pre_sur, function(x){table(x$phe$event)}))

可以看到每个癌症里面的病人的结局时间分布并不固定:

> do.call(rbind,
+         lapply(pre_sur, function(x){table(x$phe$event)}))
       0   1
ACC   51  28
BLCA 228 178
BRCA 925 151
CESC 220  71
CHOL  18  18
COAD 340  98
DLBC  38   9
ESCA  97  64

接下来进行批量单基因cox分析

必须要保证前面的表达量矩阵信息跟这个临床信息相对应哦,这样我的代码才有可能是复制粘贴就能运行:

library(survminer) 
library(survival)
cox_list = lapply( pre_sur , function(x){
  phe=x$phe
  exprSet=x$exprSet
  mySurv <- with(phe, Surv(time, event))
  
  cox_results <-apply(exprSet , 1 , function(gene){
    # gene= as.numeric(exprSet[1,])
    group=ifelse(gene>median(gene),'high','low'
    if( length(table(group))<2)
      return(NULL)
    survival_dat <- data.frame(group=group,# stage=phe$stage,
                               stringsAsFactors = F)
    m=coxph(mySurv ~ group, 
            # mySurv ~  stage+ group,  # 如果有交叉变量
            data =  survival_dat)
    
    beta <- coef(m)
    se <- sqrt(diag(vcov(m)))
    HR <- exp(beta)
    HRse <- HR * se
    
    #summary(m)
    tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^21),
                       HR = HR, HRse = HRse,
                       HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^21),
                       HRCILL = exp(beta - qnorm(.97501) * se),
                       HRCIUL = exp(beta + qnorm(.97501) * se)), 3)
    return(tmp['grouplow',])
    
  })
  return(t(cox_results))
  
})
names(cox_list)
names(cox_list)=gsub('TCGA-','',gsub('.htseq_counts.Rdata','',fs)) 
save(cox_list,file = 'batch_cox_list.Rdata')

我们简单的看看其中一个癌症的结果矩阵:

> load(file = 'batch_cox_list.Rdata')
> head(cox_list[[1]])
         coef   se     z     p   HR HRse   HRz  HRp HRCILL HRCIUL
ZZZ3   -0.470 0.38 -1.23 0.220 0.62 0.24 -1.57 0.12   0.29   1.32
ZZEF1  -0.421 0.38 -1.10 0.273 0.66 0.25 -1.36 0.17   0.31   1.39
ZYX    -0.224 0.38 -0.58 0.559 0.80 0.31 -0.66 0.51   0.38   1.70
ZYG11B -0.051 0.38 -0.13 0.894 0.95 0.36 -0.14 0.89   0.45   2.01
ZYG11A -1.235 0.40 -3.12 0.002 0.29 0.12 -6.17 0.00   0.13   0.63
ZXDC   -0.212 0.38 -0.56 0.577 0.81 0.31 -0.62 0.54   0.38   1.70

要根据这个批量单基因cox分析结果来挑选保护因子和风险因子,跟前面的使用了固定阈值,  logFC_cutoff = 1 以及  adj.P.Val<0.05,确定每个癌症的统计学显著的表达量上下调基因是一回事。

我们仍然是需要首先筛选p值确定统计学显著,然后根据HR值来判断保护因子和风险因子:

关于HR值:

  • HR = 1: No effect
  • HR < 1: Reduction in the hazard
  • HR > 1: Increase in Hazard

Note that in cancer studies:

  • A covariate with hazard ratio > 1 (i.e.: b > 0) is called bad prognostic factor
  • A covariate with hazard ratio < 1 (i.e.: b < 0) is called good prognostic factor

值得注意的是理论上HR > 1是增加了肿瘤的风险,意思是该基因表达量升高会增加肿瘤风险,但是我上面的代码里面选取的是     return(tmp['grouplow',]) ,所以我们的结果HR > 1的含义是该基因表达量降低会增加肿瘤风险。

解释起来会有一点点绕,不过这样的结果很容易通过一个分组KM曲线去肉眼检查一下基因具体到底是保护因子还是风险因子,就跟我们肉眼检查表达量上下调基因会使用箱线图一样。

跟表达量上下调基因

跟前面的笔记:癌基因一定在肿瘤部位高表达吗 的代码大同小异:

load(file = 'batch_cox_list.Rdata')
head(cox_list[[1]]) 
library(data.table)
drivers = fread('canonical_drivers.txt',data.table = F)[,1]
tsg = fread('Human_TSGs.txt',data.table = F)[,2]
length(intersect(drivers,tsg))
head(intersect(drivers,tsg))

sum1 = lapply(names(cox_list), function(x){
  # x= names(cox_list)[[11]]
  print(x)
  cox_results = as.data.frame(cox_list[[x]])
 if(nrow(cox_results)<10){
   return(NULL)
 }else{
   head(cox_results) 
   cox_results=cox_results[cox_results$p < 0.05 ,]
   gene_up= rownames( cox_results[ cox_results$HR > 1 ,])
   gene_down= rownames( cox_results[ cox_results$HR < 1 ,]) 
   
   df=data.frame(
     cancer=x,
     group=c('up','down'),
     deg=c(length(gene_up),length(gene_down)),
     inter1 = c(length(intersect(gene_up,drivers)),length(intersect(gene_down,drivers))),
     inter2 = c(length(intersect(gene_up,tsg)),length(intersect(gene_down,tsg))) 
   )
   return(df)
 }
})
sum1 = do.call(rbind,sum1 )
sum1$ratio1 = sum1$inter1/length(drivers)
sum1$ratio2 = sum1$inter2/length(tsg)
options(digits = 2)
sum1

结果还是蛮让我意外的,如下所示:

> sum1
   cancer group  deg inter1 inter2 ratio1  ratio2
1     ACC    up  606     14     37 0.0237 0.03040
2     ACC  down 4338    150    219 0.2542 0.17995
3    BLCA    up 2464     67    103 0.1136 0.08463
4    BLCA  down  696     14     60 0.0237 0.04930
5    BRCA    up  775     33     57 0.0559 0.04684
6    BRCA  down  700     16     36 0.0271 0.02958
7    CESC    up 2141     71    101 0.1203 0.08299
8    CESC  down  245      5     20 0.0085 0.01643
9    CHOL    up  427     19     21 0.0322 0.01726
10   CHOL  down  103      3      7 0.0051 0.00575
11   COAD    up   95      3      5 0.0051 0.00411
12   COAD  down  586     14     21 0.0237 0.01726
13   DLBC    up   79      1      5 0.0017 0.00411
14   DLBC  down   73      1      1 0.0017 0.00082

因为每个癌症都是使用了同样的阈值,所以各自的保护因子和风险因子数量不一样。但是可以看到,跟前面的笔记:癌基因一定在肿瘤部位高表达吗 的结论类似,并没有明显的倾向性。

肯定是有一些癌基因,它在肿瘤里面高表达却不影响生存,但是并不是全部的癌基因哦,因为导致肿瘤发生发展只需要少量关键基因即可。

同理肯定是有一些抑癌基因在肿瘤里面表达量降低,但也没有影响生存。

其实生存分析受到了的干扰因素非常多,一个目标基因可能是非常有临床意义所以它有统计学显著的生存意义,但是两万多个基因总是有那么一些基因跟目标基因表达量相关性非常高,所以也有 统计学显著的生存意义 。但是我们确实不能因为某个基因有统计学显著的生存意义就判定他一定是很重要。

比如,我们随机选取几个基因看看:

> head(cox_list[[1]][cox_list[[1]][,4]<0.05,])
        coef   se    z     p    HR  HRse   HRz   HRp HRCILL HRCIUL
ZYG11A -1.24 0.40 -3.1 0.002 0.291 0.115  -6.2 0.000  0.134   0.63
ZWINT  -2.38 0.54 -4.4 0.000 0.092 0.050 -18.1 0.000  0.032   0.27
ZWILCH -1.75 0.47 -3.8 0.000 0.174 0.081 -10.2 0.000  0.070   0.43
ZW10   -1.11 0.42 -2.6 0.008 0.329 0.138  -4.9 0.000  0.145   0.75
ZSWIM4 -0.83 0.40 -2.1 0.037 0.436 0.174  -3.2 0.001  0.200   0.95
ZSWIM3  0.92 0.41  2.3 0.024 2.505 1.016   1.5 0.139  1.131   5.55

批量绘制多个基因的KM曲线代码是:


## 针对第一个癌症,具体检查几个基因
head(cox_list[[1]][cox_list[[1]][,4]<0.05,])
cg =  head(rownames(cox_list[[1]][cox_list[[1]][,4]<0.05,]))
cg 
x=pre_sur[[1]]
phe=x$phe
exprSet=x$exprSet 
identical(colnames(exprSet), rownames(phe)) 
mySurv <- with(phe, Surv(time, event))
library(ggstatsplot)
overlap_list <- lapply(cg, function(i){
  # i = cg[1]
  survival_dat = phe
  gene = as.numeric(exprSet[i,])
  survival_dat$gene = ifelse(gene > median(gene),'high','low')
  table(survival_dat$gene)
  library(survival)
  fit <- survfit(Surv(time, event) ~ gene,
                 data = survival_dat)
  
  survp = ggsurvplot(fit,data = survival_dat, #这里很关键,不然会报错
                     legend.title = i, #定义图例的名称 
                     # legend = "top",#图例位置
                     # legend.labs = c('High', 'Low'),
                     pval = T#在图上添加log rank检验的p值
                     # pval.method = TRUE,#添加p值的检验方法
                     risk.table = TRUE
                     risk.table.y.text = F,
                     xlab = "Time in years"#x轴标题
                     # xlim = c(0, 10), #展示x轴的范围
                     # break.time.by = 1, #x轴间隔
                     size = 1.5#线条大小
                     ggtheme = theme_ggstatsplot(),
                     palette="nejm"#配色
  )
  return(survp)               
}) 

n = length(overlap_list)
n
x=floor(n/2);y=2
all_plot <- arrange_ggsurvplots(overlap_list,print = F,ncol =x, nrow = y,
                                risk.table.height = 0.3,
                                surv.plot.height = 0.7)

all_plot  

出图如下:

批量绘制多个基因的KM曲线代码

我在生信技能树多次分享过生存分析的细节;

生存分析是目前肿瘤等疾病研究领域的点睛之笔!

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存